##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(dplyr)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--cor_file"), type = "character"),
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata"

    ## 单细胞的atac
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## 表达和motif的相关性
    cor_file <- "~/20231121_singleMuti/results/celltype_plot/mfuzz/cor.motif_atac-rna.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v3/germ/"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
out_path <- opt$out_path
cor_file <- opt$cor_file
comine_data_file <- opt$comine_data_file

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
a <- load(rna_file)
b <- load(comine_data_file)
dat_cor <- data.frame(fread(cor_file))

###########################################################################################
## 表达
exp_numeric <- apply(scrnat@assays$RNA@data , 1 , as.numeric)
rownames(exp_numeric) <- colnames(scrnat@assays$RNA@data)
exp_numeric <- t(exp_numeric)
colnames(exp_numeric) <- paste0(scrnat$cell , "_" , scrnat$seurat_clusters)

out_file <- paste0(out_path , "/GeneExpression.All.rds")
saveRDS( exp_numeric , out_file )

###########################################################################################
## motif
motif_se <- getMatrixFromProject(testis_combined_peak_combineRNA, useMatrix="MotifMatrix")
motif_mat <- assays(motif_se)$deviation
## 转化为数值矩阵
motif_mat_num <- apply(motif_mat , 1 , as.numeric)
rownames(motif_mat_num) <- colnames(motif_mat)
motif_mat_num <- t(motif_mat_num)
## 去除ENSG基因
#motif_mat_num <- motif_mat_num[grep( "ENSG" , rownames(motif_mat_num) , invert = T),]
## 去除没有表达的motif
motif_mat_num <- motif_mat_num[dat_cor$MotifMatrix_name,]
## 去除名字的_
rownames(motif_mat_num) <- sapply(strsplit(rownames(motif_mat_num) , "_") , "[" , 1)

cell_cluster <- data.frame( cell = rownames(testis_combined_peak_combineRNA@cellColData) , cluster = testis_combined_peak_combineRNA@cellColData$seurat_clusters )
cell_cluster$use <- paste0( cell_cluster$cell , "_" , cell_cluster$cluster )
rownames(cell_cluster) <- cell_cluster$cell

## 用于计算motif的评价表达
GSM_mat_num <- motif_mat_num

colnames(motif_mat_num) <- cell_cluster[colnames(motif_mat_num),"use"]

out_file <- paste0(out_path , "/Motif.All.rds")
saveRDS( motif_mat_num , out_file )

## 计算motif在每个细胞的评价表达
sco_motif <- sapply(unique(unique(testis_combined_peak_combineRNA@cellColData$cell_type)),function(x){
    print(x)
    sapply(unique(rownames(GSM_mat_num)),function(y){
        mean(
            as.numeric(
                as.vector(GSM_mat_num[y,which(testis_combined_peak_combineRNA@cellColData$cell_type==x)])
            )
        )
    })
})

sco_motif <- data.frame(sco_motif)
sco_motif$gene <- rownames(sco_motif)
sco_motif <- sco_motif[,c( ncol(sco_motif) , (1:ncol(sco_motif)-1) )]

out_file <- paste0(out_path , "/Motif.MeanByCellType.tsv")
write.table( sco_motif , out_file , row.names = F , quote = F , sep = "\t" )


